/*******************************************************************************
 * NAME
 *        geom_inverse - geometry transformation inverse function
 *
 * SYNOPSIS
 *        #include <stdio.h>
 *        #include <string.h>
 *        #include <mcce.h>
 *
 *        void geom_inverse(GEOM *op);
 *
 * DESCRIPTION
 *        The geom_inverse() function reverses the transformation recorded in the
 *        operation recorder op by inverting the passed in transformation matrix.
 *
 * SEE ALSO
 *        geom_reset, geom_dump, geom_move, geom_roll, geom_3v_onto_3v
 *
 * EXAMPLE
 *        #include <stdio.h>
 *        #include <string.h>
 *        #include <mcce.h>
 *
 *        int main()
 *        {  GEOM recorder;
 *           VECTOR v;
 *
 *           v.x = 1.0; v.y = 1.0; v.z = 1.0;
 *           geom_reset(&recorder);
 *           geom_move(&recorder, v);
 *           geom_inverse(&op);
 *           geom_dump(recorder);
 *
 *           return 0;
 *        }
 *
 * DEPENDENCY
 *        det4
 *
 * AUTHOR
 *        Junjun Mao, 05/27/2003
 *******************************************************************************/

#include <stdio.h>
#include <string.h>
#include "mcce.h"

void geom_inverse(GEOM *op)
{  float t[4][4];
   float s;		/* the determinant of this matrix */

   memcpy(t, op->M, 16*sizeof(float));

   s = det4(t);

   op->M[0][0] = (t[1][1] * (t[2][2] * t[3][3] - t[2][3] * t[3][2]) +
                  t[1][2] * (t[2][3] * t[3][1] - t[2][1] * t[3][3]) +
                  t[1][3] * (t[2][1] * t[3][2] - t[2][2] * t[3][1])) / s;

   op->M[0][1] = (t[2][1] * (t[0][2] * t[3][3] - t[0][3] * t[3][2]) +
                  t[2][2] * (t[0][3] * t[3][1] - t[0][1] * t[3][3]) +
                  t[2][3] * (t[0][1] * t[3][2] - t[0][2] * t[3][1])) / s;

   op->M[0][2] = (t[3][1] * (t[0][2] * t[1][3] - t[0][3] * t[1][2]) +
                  t[3][2] * (t[0][3] * t[1][1] - t[0][1] * t[1][3]) +
                  t[3][3] * (t[0][1] * t[1][2] - t[0][2] * t[1][1])) / s;

   op->M[0][3] = (t[0][1] * (t[1][3] * t[2][2] - t[1][2] * t[2][3]) +
                  t[0][2] * (t[1][1] * t[2][3] - t[1][3] * t[2][1]) +
                  t[0][3] * (t[1][2] * t[2][1] - t[1][1] * t[2][2])) / s;

   op->M[1][0] = (t[1][2] * (t[2][0] * t[3][3] - t[2][3] * t[3][0]) +
                  t[1][3] * (t[2][2] * t[3][0] - t[2][0] * t[3][2]) +
                  t[1][0] * (t[2][3] * t[3][2] - t[2][2] * t[3][3])) / s;

   op->M[1][1] = (t[2][2] * (t[0][0] * t[3][3] - t[0][3] * t[3][0]) +
                  t[2][3] * (t[0][2] * t[3][0] - t[0][0] * t[3][2]) +
                  t[2][0] * (t[0][3] * t[3][2] - t[0][2] * t[3][3])) / s;

   op->M[1][2] = (t[3][2] * (t[0][0] * t[1][3] - t[0][3] * t[1][0]) +
                  t[3][3] * (t[0][2] * t[1][0] - t[0][0] * t[1][2]) +
                  t[3][0] * (t[0][3] * t[1][2] - t[0][2] * t[1][3])) / s;

   op->M[1][3] = (t[0][2] * (t[1][3] * t[2][0] - t[1][0] * t[2][3]) +
                  t[0][3] * (t[1][0] * t[2][2] - t[1][2] * t[2][0]) +
                  t[0][0] * (t[1][2] * t[2][3] - t[1][3] * t[2][2])) / s;

   op->M[2][0] = (t[1][3] * (t[2][0] * t[3][1] - t[2][1] * t[3][0]) +
                  t[1][0] * (t[2][1] * t[3][3] - t[2][3] * t[3][1]) +
                  t[1][1] * (t[2][3] * t[3][0] - t[2][0] * t[3][3])) / s;

   op->M[2][1] = (t[2][3] * (t[0][0] * t[3][1] - t[0][1] * t[3][0]) +
                  t[2][0] * (t[0][1] * t[3][3] - t[0][3] * t[3][1]) +
                  t[2][1] * (t[0][3] * t[3][0] - t[0][0] * t[3][3])) / s;

   op->M[2][2] = (t[3][3] * (t[0][0] * t[1][1] - t[0][1] * t[1][0]) +
                  t[3][0] * (t[0][1] * t[1][3] - t[0][3] * t[1][1]) +
                  t[3][1] * (t[0][3] * t[1][0] - t[0][0] * t[1][3])) / s;

   op->M[2][3] = (t[0][3] * (t[1][1] * t[2][0] - t[1][0] * t[2][1]) +
                  t[0][0] * (t[1][3] * t[2][1] - t[1][1] * t[2][3]) +
                  t[0][1] * (t[1][0] * t[2][3] - t[1][3] * t[2][0])) / s;

   op->M[3][0] = (t[1][0] * (t[2][2] * t[3][1] - t[2][1] * t[3][2]) +
                  t[1][1] * (t[2][0] * t[3][2] - t[2][2] * t[3][0]) +
                  t[1][2] * (t[2][1] * t[3][0] - t[2][0] * t[3][1])) / s;

   op->M[3][1] = (t[2][0] * (t[0][2] * t[3][1] - t[0][1] * t[3][2]) +
                  t[2][1] * (t[0][0] * t[3][2] - t[0][2] * t[3][0]) +
                  t[2][2] * (t[0][1] * t[3][0] - t[0][0] * t[3][1])) / s;

   op->M[3][2] = (t[3][0] * (t[0][2] * t[1][1] - t[0][1] * t[1][2]) +
                  t[3][1] * (t[0][0] * t[1][2] - t[0][2] * t[1][0]) +
                  t[3][2] * (t[0][1] * t[1][0] - t[0][0] * t[1][1])) / s;

   op->M[3][3] = (t[0][0] * (t[1][1] * t[2][2] - t[1][2] * t[2][1]) +
                  t[0][1] * (t[1][2] * t[2][0] - t[1][0] * t[2][2]) +
                  t[0][2] * (t[1][0] * t[2][1] - t[1][1] * t[2][0])) / s;
   return;
}
